JP32 (ActD CP mouse): BAM reclustering

rr suppressPackageStartupMessages(library(dplyr)) suppressPackageStartupMessages(library(Seurat)) suppressPackageStartupMessages(library(clustree)) suppressPackageStartupMessages(library(ggplot2)) suppressPackageStartupMessages(library(VennDiagram)) suppressPackageStartupMessages(library(cowplot)) suppressPackageStartupMessages(library(harmony))

rr path<-~/Documents/Kia/CPepi_review/results/
sample.names<-c(32) subset.name<-
source(~/Documents/scripts/preprosessing_May.2020.R) gene.markers<-read.csv(~/Documents/gene_protein_lists/markergenes_mouse.csv) cols.use=c(1,3,1,,,1,,,3,,4,,,,3,1,, , ,, 3, 30, , , , 30)

rr seur_full<-readRDS(paste0(path,sample.names,.seurat.no.doublets.rds))

rr DimPlot(seur_full,label=T, repel=T,cols=cols.use, group.by = ,reduction = 30)+ggtitle( dataset : annotation)

rr DimPlot(seur_full,label=T, repel=T,cols=cols.use, group.by = _cell_type,reduction = 30)+ggtitle(dataset : Manuscript cell types)

Subset BAMs

rr Idents(seur_full)=
clusters<-c(paste(, 1:5), prolif,1,2 ,  ) clusters[!clusters %in% Idents(seur_full)] #check if some of the cluster names is not correct

character(0)

Subset the cells

rr seur<-seur_full[,WhichCells(seur_full, idents=clusters)] seur$full.dataset.annot=Idents(seur)

Remove genes that are not expressed in any cell

rr num.cells.per.gene <- rowSums(as.matrix(GetAssayData(seur, slot = )) > 0) genes.use <- names(num.cells.per.gene[which(num.cells.per.gene >= 1)]) seur<-subset(seur,features=genes.use)

rr cat(/cells of full data : ,dim(seur_full))

genes/cells of full data : 
 12283 2155

rr cat(\ngenes/cells of subset : ,dim(seur))

genes/cells of subset : 
 12177 1486

Standard preprocessing and UMAP of the subsetted dataset (Selecting variable genes, scaling, PCA, UMAP )

rr seur <- NormalizeData(seur,verbose = F) seur <- FindVariableFeatures(seur,verbose=F) seur <- ScaleData(seur,verbose = F) seur <- RunPCA(seur, features = VariableFeatures(seur),verbose=F)

rr ElbowPlot(object = seur,ndims =50)

PC 1:10

rr dims.use=10

rr seur <- RunUMAP(seur, dims = 1:dims.use, verbose=F,reduction.name =paste0(,dims.use),reduction.key =paste0(,dims.use,_))

rr DimPlot(object = seur, group.by = _cell_type,label=T, repel=T,reduction=paste0(,dims.use))+ggtitle(paste(subset.name,: Manuscript cell types))

rr DimPlot(object = seur, group.by = .dataset.annot,label=T, repel=T,reduction=paste0(,dims.use))+ggtitle(paste(subset.name,: annotation from full dataset))

rr DimPlot(object = seur, group.by = .immgen.main,label=T, repel=T, cols=cols.use,reduction=paste0(,dims.use))+ggtitle(paste(subset.name,: singler.immgen.main))

rr FeaturePlot(seur, c(6c2,2,4a7,2,,,2ry12,7, 3e, 67, 5, 3),reduction=paste0(,dims.use))

Clustering with Leiden algorithm

rr seur <- FindNeighbors(seur, dims = 1:dims.use, verbose=F,graph.name =paste0(_snn_PC,dims.use)) for ( i in seq(0,2, 0.25)) seur <- FindClusters(seur, resolution = i, verbose=F, algorithm = 4,graph.name =paste0(_snn_PC,dims.use)) # algorithm= 4 is Leiden algorithm - often performs better

Plot of a clustering tree showing the relationship between clusterings at different resolutions. (using the clustree package)

rr clustree(seur, prefix = paste0(_snn_PC,dims.use,_res.)) + ggtitle(paste(subset.name,: Clustering tree PC=, dims.use))

rr plot<-list() for ( res in c(0.5, 0.75,1,1.25)) plot[[as.character(res)]]<-DimPlot(seur, pt.size = 1,label=T,repel=T, group.by = paste0(_snn_PC,dims.use,_res.,res),reduction=paste0(,dims.use)) + ggtitle(paste(=,dims.use,=,res)) plot_grid(plotlist=plot)

Let’s find differentially expressed genes per cluster

rr res=0.75 Idents(seur)= paste0(_snn_PC,dims.use,_res.,res)

rr res=0.75 dims.use=10 DEgenes_list <- readRDS( paste0( path,.res,res,_,,dims.use,_,subset.name,_,sample.names,.rds))

rr features.use=unlist(lapply(DEgenes_list, function(x) { head(x[x$avg_logFC>0,]$gene)})) DoHeatmap(seur, features = features.use, assay = , angle = 90, label =T, size=4) + scale_fill_gradient2(low = , mid = ,high = )+ theme(axis.text.y= element_text(size=11))+ ggtitle(paste( res=,res,  PC= , dims.use))

Scale for 'fill' is already present. Adding another scale for 'fill', which will replace the existing scale.

rr DimPlot(object = seur,label=T, repel=T, reduction=paste0(,dims.use))+ggtitle(paste( res=,res,  PC= , dims.use))

rr ### Save as excel table setwd(path) first_sheet_name<-names(DEgenes_list)[1] #JAVA specific garbage collection jgc <- function() { rJava::.jcall(/lang/System, method = ) } #Create the excel file and add the first sheet write.xlsx2(DEgenes_list[[first_sheet_name]], file=paste0(.res,res,_,,dims.use,_,subset.name,_,sample.names,.xlsx), sheetName=first_sheet_name, row.names=FALSE) #Add the remaining sheets tot he excel file for ( i in names(DEgenes_list)[names(DEgenes_list)!=first_sheet_name]) { gc() jgc() message(sheet , i) write.xlsx2(DEgenes_list[[i]], file=paste0(.res,res,_,,dims.use,_,subset.name,_,sample.names,.xlsx), sheetName=i, append=TRUE, row.names=FALSE) }

---

PC 1:15

rr dims.use=15

rr seur <- RunUMAP(seur, dims = 1:dims.use, verbose=F,reduction.name =paste0(,dims.use),reduction.key =paste0(,dims.use,_))

rr DimPlot(object = seur, group.by = _cell_type,label=T, repel=T,reduction=paste0(,dims.use))+ggtitle(paste(subset.name,: Manuscript cell types))

rr DimPlot(object = seur, group.by = .dataset.annot,label=T, repel=T,reduction=paste0(,dims.use))+ggtitle(paste(subset.name,: annotation from full dataset))

Clustering with Leiden algorithm

rr seur <- FindNeighbors(seur, dims = 1:dims.use, verbose=F,graph.name =paste0(_snn_PC,dims.use)) for ( i in seq(0,2, 0.25)) seur <- FindClusters(seur, resolution = i, verbose=F, algorithm = 4,graph.name =paste0(_snn_PC,dims.use)) # algorithm= 4 is Leiden algorithm - often performs better

Plot of a clustering tree showing the relationship between clusterings at different resolutions. (using the clustree package)

rr clustree(seur, prefix = paste0(_snn_PC,dims.use,_res.)) + ggtitle(paste(subset.name,: Clustering tree PC=, dims.use))

rr plot<-list() for ( res in c(0.5, 0.75,1,1.25)) plot[[as.character(res)]]<-DimPlot(seur, pt.size = 1,label=T,repel=T, group.by = paste0(_snn_PC,dims.use,_res.,res),reduction=paste0(,dims.use)) + ggtitle(paste(=,dims.use,=,res)) plot_grid(plotlist=plot)

Visualize the same clusterings on UMAP with 10 PC

rr plot<-list() for ( res in c(0.5, 0.75,1,1.25)) plot[[as.character(res)]]<-DimPlot(seur,label=T,repel=T, group.by = paste0(_snn_PC,dims.use,_res.,res), reduction= paste0(10)) + ggtitle(paste(=,dims.use,=,res)) plot_grid(plotlist=plot)

Let’s find differentially expressed genes per cluster

rr res=0.75 Idents(seur)= paste0(_snn_PC,dims.use,_res.,res)

rr res=0.75 dims.use=15 DEgenes_list <- readRDS( paste0( path,.res,res,_,,dims.use,_,subset.name,_,sample.names,.rds))

rr features.use=unlist(lapply(DEgenes_list, function(x) { head(x[x$avg_logFC>0,]$gene)})) DoHeatmap(seur, features = features.use, assay = , angle = 90, label =T, size=4) + scale_fill_gradient2(low = , mid = ,high = )+ theme(axis.text.y= element_text(size=11))+ ggtitle(paste( res=,res,  PC= , dims.use))

Scale for 'fill' is already present. Adding another scale for 'fill', which will replace the existing scale.

rr DimPlot(object = seur,label=T, repel=T, reduction=paste0(,dims.use))+ggtitle(paste( res=,res,  PC= , dims.use))

rr ### Save as excel table setwd(path) first_sheet_name<-names(DEgenes_list)[1] #JAVA specific garbage collection jgc <- function() { rJava::.jcall(/lang/System, method = ) } #Create the excel file and add the first sheet write.xlsx2(DEgenes_list[[first_sheet_name]], file=paste0(.res,res,_,,dims.use,_,subset.name,_,sample.names,.xlsx), sheetName=first_sheet_name, row.names=FALSE) #Add the remaining sheets tot he excel file for ( i in names(DEgenes_list)[names(DEgenes_list)!=first_sheet_name]) { gc() jgc() message(sheet , i) write.xlsx2(DEgenes_list[[i]], file=paste0(.res,res,_,,dims.use,_,subset.name,_,sample.names,.xlsx), sheetName=i, append=TRUE, row.names=FALSE) }

---

PC 1:20

rr dims.use=20

rr seur <- RunUMAP(seur, dims = 1:dims.use, verbose=F,reduction.name =paste0(,dims.use),reduction.key =paste0(,dims.use,_))

rr DimPlot(object = seur, group.by = _cell_type,label=T, repel=T,reduction=paste0(,dims.use))+ggtitle(paste(subset.name,: Manuscript cell types))

rr DimPlot(object = seur, group.by = .dataset.annot,label=T, repel=T,reduction=paste0(,dims.use))+ggtitle(paste(subset.name,: annotation from full dataset))

Clustering with Leiden algorithm

rr seur <- FindNeighbors(seur, dims = 1:dims.use, verbose=F,graph.name =paste0(_snn_PC,dims.use)) for ( i in seq(0,2, 0.25)) seur <- FindClusters(seur, resolution = i, verbose=F, algorithm = 4,graph.name =paste0(_snn_PC,dims.use)) # algorithm= 4 is Leiden algorithm - often performs better

Plot of a clustering tree showing the relationship between clusterings at different resolutions. (using the clustree package)

rr clustree(seur, prefix = paste0(_snn_PC,dims.use,_res.)) + ggtitle(paste(subset.name,: Clustering tree PC=, dims.use))

rr plot<-list() for ( res in c(0.5, 0.75,1,1.25)) plot[[as.character(res)]]<-DimPlot(seur, pt.size = 1,label=T,repel=T, group.by = paste0(_snn_PC,dims.use,_res.,res),reduction=paste0(,dims.use)) + ggtitle(paste(=,dims.use,=,res)) plot_grid(plotlist=plot)

Visualize the same clusterings on UMAP with 10 PC

rr plot<-list() for ( res in c(0.5, 0.75,1,1.25)) plot[[as.character(res)]]<-DimPlot(seur,label=T,repel=T, group.by = paste0(_snn_PC,dims.use,_res.,res), reduction= paste0(10)) + ggtitle(paste(=,dims.use,=,res)) plot_grid(plotlist=plot)

Let’s find differentially expressed genes per cluster

rr res=1 Idents(seur)= paste0(_snn_PC,dims.use,_res.,res)

rr res=1 dims.use=10 DEgenes_list <- readRDS( paste0( path,.res,res,_,,dims.use,_,subset.name,_,sample.names,.rds))

rr features.use=unlist(lapply(DEgenes_list, function(x) { head(x[x$avg_logFC>0,]$gene)})) DoHeatmap(seur, features = features.use, assay = , angle = 90, label =T, size=4) + scale_fill_gradient2(low = , mid = ,high = )+ theme(axis.text.y= element_text(size=10))+ ggtitle(paste( res=,res,  PC= , dims.use))

Scale for 'fill' is already present. Adding another scale for 'fill', which will replace the existing scale.

rr ### Save as excel table setwd(path) first_sheet_name<-names(DEgenes_list)[1] #JAVA specific garbage collection jgc <- function() { rJava::.jcall(/lang/System, method = ) } #Create the excel file and add the first sheet write.xlsx2(DEgenes_list[[first_sheet_name]], file=paste0(.res,res,_,,dims.use,_,subset.name,_,sample.names,.xlsx), sheetName=first_sheet_name, row.names=FALSE) #Add the remaining sheets tot he excel file for ( i in names(DEgenes_list)[names(DEgenes_list)!=first_sheet_name]) { gc() jgc() message(sheet , i) write.xlsx2(DEgenes_list[[i]], file=paste0(.res,res,_,,dims.use,_,subset.name,_,sample.names,.xlsx), sheetName=i, append=TRUE, row.names=FALSE) }

rr saveRDS(seur,paste0(path,subset.name,sample.names,.seurat.rds))

rr seur<-readRDS(paste0(path,subset.name,sample.names,.seurat.rds))

rr sessionInfo()

R version 4.0.3 (2020-10-10)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 20.04.2 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0
LAPACK: /home/daliya/anaconda3/lib/libmkl_rt.so

locale:
 [1] LC_CTYPE=en_US.UTF-8          LC_NUMERIC=C                  LC_TIME=de_BE.UTF-8           LC_COLLATE=en_US.UTF-8       
 [5] LC_MONETARY=de_BE.UTF-8       LC_MESSAGES=en_US.UTF-8       LC_PAPER=de_BE.UTF-8          LC_NAME=de_BE.UTF-8          
 [9] LC_ADDRESS=de_BE.UTF-8        LC_TELEPHONE=de_BE.UTF-8      LC_MEASUREMENT=de_BE.UTF-8    LC_IDENTIFICATION=de_BE.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] enrichR_2.1    biomaRt_2.46.0 xlsx_0.6.5     harmony_1.0    Rcpp_1.0.6     shiny_1.5.0    cowplot_1.1.1  ggrepel_0.9.1  clustree_0.4.3
[10] ggraph_2.0.4   ggplot2_3.3.3  Seurat_3.2.3   dplyr_1.0.3   

loaded via a namespace (and not attached):
  [1] reticulate_1.18             R.utils_2.10.1              tidyselect_1.1.0            RSQLite_2.2.2               AnnotationDbi_1.52.0       
  [6] htmlwidgets_1.5.3           grid_4.0.3                  BiocParallel_1.24.1         Rtsne_0.15                  devtools_2.3.2             
 [11] DropletUtils_1.10.2         munsell_0.5.0               codetools_0.2-16            ica_1.0-2                   future_1.21.0              
 [16] miniUI_0.1.1.1              withr_2.4.0                 colorspace_2.0-0            Biobase_2.50.0              knitr_1.30                 
 [21] rstudioapi_0.13             stats4_4.0.3                SingleCellExperiment_1.12.0 ROCR_1.0-11                 tensor_1.5                 
 [26] rJava_0.9-13                listenv_0.8.0               MatrixGenerics_1.2.0        labeling_0.4.2              GenomeInfoDbData_1.2.4     
 [31] polyclip_1.10-0             bit64_4.0.5                 farver_2.0.3                rhdf5_2.34.0                rprojroot_2.0.2            
 [36] parallelly_1.23.0           vctrs_0.3.6                 generics_0.1.0              xfun_0.20                   BiocFileCache_1.14.0       
 [41] R6_2.5.0                    GenomeInfoDb_1.26.2         graphlayouts_0.7.1          rsvd_1.0.3                  locfit_1.5-9.4             
 [46] AnnotationFilter_1.14.0     bitops_1.0-6                rhdf5filters_1.2.0          spatstat.utils_2.1-0        DelayedArray_0.16.0        
 [51] assertthat_0.2.1            promises_1.1.1              scales_1.1.1                gtable_0.3.0                beachmat_2.6.4             
 [56] Cairo_1.5-12.2              globals_0.14.0              processx_3.4.5              goftest_1.2-2               ensembldb_2.14.0           
 [61] tidygraph_1.2.0             rlang_0.4.10                splines_4.0.3               rtracklayer_1.50.0          lazyeval_0.2.2             
 [66] checkmate_2.0.0             yaml_2.2.1                  reshape2_1.4.4              abind_1.4-5                 backports_1.2.1            
 [71] GenomicFeatures_1.42.1      httpuv_1.5.5                usethis_2.0.0               tools_4.0.3                 ellipsis_0.3.1             
 [76] RColorBrewer_1.1-2          BiocGenerics_0.36.0         sessioninfo_1.1.1           ggridges_0.5.3              plyr_1.8.6                 
 [81] sparseMatrixStats_1.2.0     progress_1.2.2              zlibbioc_1.36.0             purrr_0.3.4                 RCurl_1.98-1.2             
 [86] ps_1.5.0                    prettyunits_1.1.1           rpart_4.1-15                openssl_1.4.3               deldir_0.2-9               
 [91] pbapply_1.4-3               viridis_0.5.1               S4Vectors_0.28.1            zoo_1.8-8                   SummarizedExperiment_1.20.0
 [96] cluster_2.1.0               fs_1.5.0                    magrittr_2.0.1              data.table_1.13.6           RSpectra_0.16-0            
[101] scattermore_0.7             lmtest_0.9-38               RANN_2.6.1                  ProtGenerics_1.22.0         fitdistrplus_1.1-3         
[106] matrixStats_0.58.0          pkgload_1.1.0               hms_1.0.0                   xlsxjars_0.6.1              patchwork_1.1.1            
[111] mime_0.9                    evaluate_0.14               xtable_1.8-4                XML_3.99-0.5                readxl_1.3.1               
[116] IRanges_2.24.1              gridExtra_2.3               testthat_3.0.1              compiler_4.0.3              tibble_3.0.5               
[121] KernSmooth_2.23-17          crayon_1.4.0                R.oo_1.24.0                 htmltools_0.5.1             mgcv_1.8-33                
[126] later_1.1.0.1               tidyr_1.1.2                 DBI_1.1.1                   tweenr_1.0.1                dbplyr_2.0.0               
[131] MASS_7.3-53                 rappdirs_0.3.1              Matrix_1.3-2                cli_2.2.0                   R.methodsS3_1.8.1          
[136] parallel_4.0.3              igraph_1.2.6                GenomicRanges_1.42.0        pkgconfig_2.0.3             GenomicAlignments_1.26.0   
[141] scuttle_1.0.4               plotly_4.9.3                xml2_1.3.2                  dqrng_0.2.1                 XVector_0.30.0             
[146] callr_3.5.1                 stringr_1.4.0               digest_0.6.27               sctransform_0.3.2           RcppAnnoy_0.0.18           
[151] Biostrings_2.58.0           spatstat.data_2.1-0         rmarkdown_2.6               cellranger_1.1.0            leiden_0.3.6               
[156] uwot_0.1.10                 edgeR_3.32.1                DelayedMatrixStats_1.12.2   curl_4.3                    Rsamtools_2.6.0            
[161] rjson_0.2.20                lifecycle_0.2.0             nlme_3.1-149                jsonlite_1.7.2              Rhdf5lib_1.12.0            
[166] desc_1.2.0                  fansi_0.4.2                 viridisLite_0.3.0           askpass_1.1                 limma_3.46.0               
[171] pillar_1.4.7                lattice_0.20-41             pkgbuild_1.2.0              fastmap_1.0.1               httr_1.4.2                 
[176] survival_3.2-7              remotes_2.2.0               glue_1.4.2                  spatstat_1.64-1             png_0.1-7                  
[181] bit_4.0.4                   ggforce_0.3.2               stringi_1.5.3               HDF5Array_1.18.0            blob_1.2.1                 
[186] memoise_1.1.0               irlba_2.3.3                 future.apply_1.7.0